A Transportation $L^p$ Distance for Signal Analysis

Abstract

Transport-based distances, such as the Wasserstein distance and earth mover’s distance, have been shown to be an effective tool in signal and image analysis. The success of transport-based distances is in part due to their Lagrangian nature which allows it to capture the important variations in many signal classes. However, these distances require the signal to be non-negative and normalised. Furthermore, the signals are considered as measures and compared by redistributing (transporting) them, which does not directly take into account the signal intensity. Here, we study a transport-based distance, called the $TL^p$ distance, that combines Lagrangian and intensity modelling and is directly applicable to general, non-positive and multichannelled signals. The distance can be computed by existing numerical methods. We give an overview of the basic properties of this distance and applications to classification, with multichannelled non-positive one-dimensional signals and two-dimensional images, and colour transfer.

Appendices

Appendix 1: Performance of $TL^p_\lambda$ in Classification Problems with Simple and Oscillatory Signals

We compare the performance of $TL_\lambda ^2$, $L^2$ and OT distances with respect to classification/clustering for the three classes $\{{\mathcal {C}}_i\}_{i=1,2,3}$ of signals defined in Fig. 4. We test how each distance performs by finding the smallest number of data points such that the classes ${\mathcal {C}}_i^N=\{f_i\}_{i=1}^N\subset {\mathcal {C}}_i$ are separable. For sufficiently large N the approximation $d_{H,\rho }({\mathcal {C}}_i^N,{\mathcal {C}}_j^N) \approx d_{H,\rho }({\mathcal {C}}_i,{\mathcal {C}}_j)$ is used to simplify the computation. Similarly, as a proxy for ${\mathbb {E}}R_\rho ({\mathcal {C}}_i^N)$ we use $R_\rho (\hat{{\mathcal {C}}}^N)$ where

$$\begin{aligned} \hat{{\mathcal {C}}}_i^N = \Bigg \{ f_\ell \, : \,&\ell = \ell _{\min }^i + \frac{n-1}{N-1} \left( \ell _{\min }^i-\ell _{\max }^i \right) , \\&\quad n\in \{1,2,\dots , N\} \Bigg \} \end{aligned}$$

is the uniform sample from class ${\mathcal {C}}_i$ (recall that class ${\mathcal {C}}_i$ is parameterized by $\ell \in [\ell _{\min }^i,\ell _{\max }^i]$ and with an abuse of notation we use the subscript of $f_\ell$ to denote the dependence of $\ell$).

It follows that the class separation distances and class coverage radius are approximated by

$$\begin{aligned} d_{H,L^2}^2({\mathcal {C}}_1^N,{\mathcal {C}}_2^N)&\approx \frac{\alpha }{2}&R^2_{L^2}({\mathcal {C}}_1^N)&\approx \frac{2}{N} \\ d_{H,L^2}^2({\mathcal {C}}_1^N,{\mathcal {C}}_3^N)&\approx \frac{3\alpha }{4}&R^2_{L^2}({\mathcal {C}}_2^N)&\approx \frac{1}{N} \\ d_{H,L^2}^2({\mathcal {C}}_2^N,{\mathcal {C}}_3^N)&\approx \frac{\alpha }{4}&R^2_{L^2}({\mathcal {C}}_3^N)&\approx \frac{2\alpha }{N\gamma } \\ d_{H,\mathrm {OT}}^2({\mathcal {C}}_1^N,{\mathcal {C}}_2^N)&\approx \frac{\beta ^2\alpha }{4}&R^2_{\mathrm {OT}}({\mathcal {C}}_1^N)&\approx \frac{\alpha }{N^2} \\ d_{H,\mathrm {OT}}^2({\mathcal {C}}_1^N,{\mathcal {C}}_3^N)&\approx \frac{\beta ^2\alpha }{4}&R^2_{\mathrm {OT}}({\mathcal {C}}_2^N)&\approx \frac{\alpha }{N^2} \\ d_{H,\mathrm {OT}}^2({\mathcal {C}}_2^N,{\mathcal {C}}_3^N)&\approx \frac{\alpha \gamma ^2}{8}&R^2_{\mathrm {OT}}({\mathcal {C}}_3^N)&\approx \frac{\alpha }{N^2} \\ d_{H,TL^2_\lambda }^2({\mathcal {C}}_1^N,{\mathcal {C}}_2^N)&\approx \frac{\alpha }{2}&R^2_{TL^2_\lambda }({\mathcal {C}}_1^N)&\approx \frac{\alpha ^2}{N} \end{aligned}$$

$$\begin{aligned} d_{H,TL^2_\lambda }^2({\mathcal {C}}_1^N,{\mathcal {C}}_3^N)&\approx \frac{3\alpha }{4}&R^2_{TL^2_\lambda }({\mathcal {C}}_2^N)&\approx \frac{4\alpha ^2}{N} \\ d_{H,TL^2_\lambda }^2({\mathcal {C}}_2^N,{\mathcal {C}}_3^N)&\approx \frac{\alpha }{4}&R^2_{TL^2_\lambda }({\mathcal {C}}_3^N)&\approx \frac{\alpha ^2}{N}. \end{aligned}$$

We have

$$\begin{aligned} \kappa _{12}^2(L^2;N)&\approx \frac{\alpha N}{4},&\kappa ^2_{13}(L^2;N)&\approx \frac{3\gamma N}{8}, \\ \kappa ^2_{12}(\mathrm {OT};N)&\approx \frac{\beta ^2 N}{4},&\kappa ^2_{13}(\mathrm {OT};N)&\approx \frac{\beta ^2 N^2}{4}, \\ \kappa ^2_{12}(TL_\lambda ^2;N)&\approx \frac{N}{8\alpha },&\kappa ^2_{13}(TL_\lambda ^2;N)&\approx \frac{3N}{4\alpha }, \\ \kappa ^2_{23}(L^2;N)&\approx \frac{\gamma N}{8}, \\ \kappa ^2_{23}(\mathrm {OT};N)&\approx \frac{\gamma ^2 N^2}{8}, \end{aligned}$$

$$\begin{aligned} \kappa ^2_{23}(TL_\lambda ^2;N)&\approx \frac{N}{16 \alpha }. \end{aligned}$$

Finally we can compute $N^*$,

$$\begin{aligned} N^*_{12}(L^2)&\approx \frac{4}{\alpha },&N^*_{13}(L^2)&\approx \frac{8}{3\gamma },&N^*_{23}(L^2)&\approx \frac{8}{\gamma } \\ N^*_{12}(\mathrm {OT})&\approx \frac{2}{\beta },&N^*_{13}(\mathrm {OT})&\approx \frac{2}{\beta },&N^*_{23}(\mathrm {OT})&\approx \frac{\sqrt{8}}{\gamma } \\ N^*_{12}(TL^2)&\approx \frac{\alpha }{8},&N^*_{13}(TL^2)&\approx \frac{4\alpha }{3},&N^*_{23}(TL^2)&\approx 16 \alpha \end{aligned}$$

which for $\beta >\frac{\alpha }{2}$, $\beta >\frac{3\gamma }{4}$ and $\gamma <\frac{\sqrt{2}\alpha }{8}$ implies the ordering given Sect. 4.1.

Appendix 2: Numerical Methods

In principle any numerical method for computing OT distances capable of dealing with an arbitrary cost function can be adapted to compute the $TL^p_\lambda$ distance. Here we describe two numerical methods we used in Sect. 4.

1.1 Iterative Linear Programming

Here we describe the iterative linear programming method of Oberman and Ruan [53] which we abbreviate OR. Although this method is not guaranteed to find the minimum in (3) we find it works well in practice and is easier to implement than, for example, methods due to Schmitzer [71] that provably minimise (3) but require a more advanced refinement procedure. See also [46] and references therein for a multiscale descent approach.

The linear programming problem restricted to a subset ${\mathcal {M}}\subseteq \Omega _h^2$ is

where $c_\lambda$ is given by (4). When ${\mathcal {M}}=\Omega _h^2$ then the $TL_\lambda ^p$ distance between $(f_h,\mu _h)$ and $(g_h,\nu _h)$ is the minimum to the above linear programme. Furthermore if $\pi _h$ is the minimiser in the $TL_\lambda ^p$ distance then it is also the solution to the linear programme in ($\hbox {LP}_{h}$) for any ${\mathcal {M}}$ containing the support of $\pi _h$. That is if one already knows (or can reasonably estimate) the set of nodes ${\mathcal {M}}$ for which the optimal plan is nonzero then one need only consider the linear programme on ${\mathcal {M}}$. This is advantageous when ${\mathcal {M}}$ is a much smaller set. Motivated by Proposition 3.5 we expect to be able to write the optimal plan as a map. This implies whilst $\pi _h$ has $n^2$ unknowns we only expect n of them to be nonzero.

The method proposed by OR is given in Algorithm 1. An initial discretisation scale $h_0$ is given and an estimate $\pi _{h_0}$ found for the linear programme ($\hbox {LP}_{h}$) with ${\mathcal {M}} = \Omega _{h_0}^2$. One then iteratively finds ${\mathcal {M}}_r\subseteq \Omega _{h_r}^2$, where $h_r=\frac{h_{r-1}}{2}$, to be the set of nodes defined by the following refinement procedure. Find the set of nodes for which $\pi _{h_{r-1}}$ is nonzero, add the neighbouring nodes and then project onto the refined grid $\Omega _{h_r}^2$. The optimal plan $\pi _{h_r}$ on $\Omega _{h_r}^2$ is then estimated by solving the linear programme ($\hbox {LP}_{h}$) with ${\mathcal {M}}={\mathcal {M}}_r$.

The grid $\Omega _{h_r}$ will scale as $(2^{rd}h_0^{-1})^2$. If the linear programme is run N times then at the $r^\text {th}$ step the linear programme has on the order of $2^{rd}h_0^{-1}$ variables. In particular on the last (and most expensive) step the number of variables is $O(2^{Nd}h_0^{-1})$. This compares to size $(2^{Nd}h_0^{-1})^2$ if the linear programme was run on the final grid without this refinement procedure.

1.2 Entropic Regularisation

Cuturi, in the context of computing OT distances, proposed regularizing the minimisation in (3) with entropy [16]. This was further developed by Benamou, Carlier, Cuturi, Nenna and Peyré [5], abbreviated to BCCNP, which is the method we describe here. Instead of considering the distance $TL^p_\lambda$ we consider

$$\begin{aligned} S_\epsilon = \inf _{\pi \in \Pi (\mu ,\nu )} \left\{ \sum _{i=1}^n \sum _{j=1}^n c_\lambda (x_i,x_j;f,g) \pi _{ij} - \epsilon H(\pi ) \right\} \end{aligned}$$

where $H(\pi )= - \sum _{i=1}^n \sum _{j=1}^n \pi _{ij} \log \pi _{ij}$ is the entropy. In the OT case the distance $S_\epsilon$ is also known as the Sinkhorn distance. It is a short calculation to show

$$\begin{aligned} S_\epsilon = \epsilon \inf _{\pi \in \Pi (\mu ,\nu )} \left\{ \mathrm {KL}(\pi |{\mathcal {K}}) \right\} \end{aligned}$$

where ${\mathcal {K}}_{ij} = \exp \left( -\frac{c_\lambda (x_i,x_j:f,g)}{\epsilon }\right)$ (the exponential is taken pointwise) and $\mathrm {KL}$ is the Kullback-Leibler divergence defined by

$$\begin{aligned} \mathrm {KL}(\pi |{\mathcal {K}}) = \sum _{i=1}^n \sum _{j=1}^n \pi _{ij} \log \left( \frac{\pi _{ij}}{{\mathcal {K}}_{ij}} \right) . \end{aligned}$$

It can be shown that the optimal choice of $\pi$ for $S_\epsilon$ can be written in the form $\pi ^* = \mathrm {diag}(u) {\mathcal {K}} \mathrm {diag}(v)$ where $u,v\in {\mathbb {R}}^n$ are limits, as $r \rightarrow \infty$, of the sequence

$$\begin{aligned} v^{(0)} = {\mathbb {I}}, \quad u^{(r)} = \frac{\underline{p}}{{\mathcal {K}} v^{(r)}}, \quad v^{(r+1)} = \frac{\underline{q}}{{\mathcal {K}}^\top u^{(r)}} \end{aligned}$$

and ${\underline{p}} = (p_1,\dots ,p_n)$, ${\underline{q}}=(q_1,\dots ,q_n)$ (multiplication is the usual matrix-vector multiplication, division is pointwise and $\top$ denotes the matrix transpose). The algorithm given in 2 is a special case of iterative Bregman projections and also known as the Sinkhorn algorithm.

The stopping condition proposed in [16] is to let $\pi ^{(r)} = \mathrm {diag}(u^{(r)}) {\mathcal {K}} \mathrm {diag}(v^{(r)})$ then stop when

$$\begin{aligned} \left| \frac{\sum _{i,j=1}^n {\mathcal {K}}_{ij} \pi ^{(r)}_{ij} - \epsilon H(\pi ^{(r)})}{\sum _{i,j=1}^n {\mathcal {K}}_{ij} \pi ^{(r-1)}_{ij} - \epsilon H(\pi ^{(r-1)})} - 1 \right| < 10^{-4}. \end{aligned}$$

Note that although as $\epsilon \rightarrow 0$ we will recover the unregularised $TL_\lambda ^p$ distance we also suffer numerical instability as ${\mathcal {K}}\rightarrow 0$ exponentially in $\epsilon$. These instabilities have been addressed in, for example, [14, 72].

For OT with quadratic cost $c(x,y) = |x-y|^2$ the Sinkhorn algorithm can be more efficiently implemented using Gaussian convolutions [74]. The two numerical methods described so far use the formulation of $TL^2_\lambda$ given by (3, 4) which interprets the $TL^2_\lambda$ as an OT distance between measures $\mu$ and $\nu$ for a (non-quadratic) cost function $c_\lambda (\cdot ,\cdot ;f,g)$, hence one cannot make use of previous OT methods such as [74].

However, we also recall that we can define the $TL^2_\lambda$ distance as the OT distance between measures $(f\times \mathrm {Id})_{\#} \mu$ and $(g\times \mathrm {Id})_{\#} \nu$, see (5), in which case the entropy regularized approach can be implemented using Gaussian convolutions in dimension $d+m$ (when $p=2$), where $f:\Omega \subseteq {\mathbb {R}}^d\rightarrow {\mathbb {R}}^m$. Although this means that the numerical method is based in a higher dimension we note the success of the bilateral grid method for bilateral filters that are also based on computing a Gaussian filter in a higher dimension [10, 57]. For colour images, where $m=3$ this approach may not be efficient however for $m=1$, these ideas have the potential for an improved algorithm.